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Abstract 
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or Kohn-Sham density functional theory. The canonical density matrix perturbation theory can 
be used to calculate temperature dependent response properties from the coupled perturbed self- 
consistent field equations as in density functional perturbation theory. The method is well suited to 
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cost as a function of system size for sufficiently large non-metallic materials and metals at high 
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I. INTRODUCTION 


Materials properties such as electric conductivity, magnetic susceptibility or electrical 
polarizabilities, are defined from their response to perturbations that are governed by the 
quantum nature of the electrons. The calculation of such quantum response properties rep¬ 
resents a major challenge because of the high cost involved. In traditional calculations the 
computational complexity scales cubically, 0{N^), or worse, with the number of atoms N, 
even when effective mean field models or density functional theory are used mi- By using 
the locality of the electronic solutions it is possible to reduce the computational cost for 
sufficiently large, non-metallic, materials to scale only linearly, 0{N), with the system size 
[3Hin]. Initially, the development of linear scaling electronic structure theory was aimed 
at calculating ground state properties and not until recently has the focus shifted towards 
the computationally more demanding task of calculating the quantum response. A num¬ 
ber of approaches to a quantum perturbation theory with reduced complexity have now 
been proposed and analyzed [TTH2^ . Linear scaling quantum perturbation theory has so 
far mainly concerned properties at zero electronic temperature. Here we extend the idea 
behind linear scaling density matrix perturbation theory [T6HI9] to calculations of static 
response properties valid also at finite electronic temperatures with fractional occupation 
of the states. Our proposed canonical density matrix perturbation theory, which is appli¬ 
cable within effective single-particle formulations, such as tight-binding, Hartree-Fock or 
Kohn-Sham density functional theory, can be applied to calculate temperature dependent 
response properties from the solution of the coupled perturbed self-consistent field equations 
[111231 ED as in density functional perturbation theory [21125]. The canonical density matrix 
perturbation scheme should be directly applicable in a number of existing program packages 
for linear scaling electronic structure calculations, including CONQUEST [9l[26ll27|, CP2K 
[2S|, ergo [2S1 EZ], FEMTECK [221 ED, FreeON [22], HONPAS [22], LATTE |2E1 122], 
ONETEP [22], OPEN-MX [2H], and SIESTA [2D]. While originally motivated by its ability 
to achieve linear scaling complexity, our canonical density matrix perturbation theory is 
quite general and straightforward to use with high efficiency also for material systems that 
are too small to reach the linear scaling regime. The computational kernel of the algorithm 
is centered around generalized matrix-matrix multiplications that are well known to pro¬ 
vide close to peak performance on many computer platforms using dense algebra, including 
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graphics processing units (GPU’s) [ill H2] . 

The paper is outlined as follows; hrst we present the canonical density matrix perturbation 
theory. Thereafter we show how it can be used to calculate temperature dependent free 
energy response properties, such as static polarizabilities and hyperpolarizabilities. We 
discuss the alternative of using hnite difference schemes and its potential problems. We 
conclude by discussing the capability of the canonical density matrix perturbation theory 
to reach linear scaling complexity in the computational cost. 


II. CANONICAL DENSITY MATRIX PERTURBATION THEORY 


In our density matrix perturbation theory we will use the single-particle density matrix 
and its derivatives to represent the electronic structure and its response to perturbations. 
With the density matrix formulation it is easy to utilize matrix sparsity from electronic 
nearsightedness pumaiiii] and it allows direct calculations of observables. The effective 
single-particle density matrix, P, at the electronic temperature Tg, can be calculated from 
the Hamiltonian, H, using a recursive Fermi operator expansion [^Hl - HH] . 


( 1 ) 




where the inverse temperature (3 = l/(fcBTe), /i is the chemical potential, and I is the identity 
matrix (see Appendix). Both H and P are here assumed to be matrix representations 
in an orthogonal basis. The expansion can be calculated through intermediate matrices 
Xn = XniXn-i) for 77, = 1, 2, 3,..., M, where 


Xo = MH) = - pJ), 


Xn = Pn{Xn-l) = 






( 2 ) 


In the canonical (NVT) ensemble, the chemical potential p is chosen such that the density 
matrix has the correct occupation, Tr[P] = Wcc, where N^cc is the number of occupied 
states. The recursion scheme above provides a very efficient and rapidly converging ex¬ 
pansion and the number of recursion steps M can be kept low (M < 20). Because of the 
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particular form of the Fade polynomial each iteration involves a solution of a 

system of linear equations, which is well tailored for the linear conjugate gradient method 
[151 uni SB] • The recursive expansion avoids the calculation of individual eigenvalues and 
eigenfunctions and is therefore well suited to reach linear scaling complexity in the computa¬ 
tional cost for sufficiently large non-metallic problems, which can utilize thresholded sparse 
matrix algebra |6]. 

A canonical density matrix response expansion, 

F(A) = + APT) + + ... , (3) 

where Tr[P^^'^] = 0 for k > 0, with respect to a perturbation in the Hamiltonian, 

H{X) = + ..., (4) 


can be constructed at hnite electronic temperatures, Tg > 0, based on the recursive Fermi 
operator expansion in Eqs. Q and ([^ above. The technique is given by a free energy 
generalization of the zero temperature linear scaling density matrix perturbation theory 
[islini. The idea is to transfer the perturbations up to some specihc order in each iteration 
step in the recursive Fermi-operator expansion, i.e. 
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+ Axr_', +...) 


( 1 ) 


A=0 


(5) 


for n = 0,1,..., M, where The additional problem of conserving the number 

of particles in a canonical ensemble, which requires Tr[P^^'>] = 0 for fc > 0, is achieved by 
including the corresponding perturbative expansion of the chemical potential, i.e. 


yU = /i(A) = -I- ApT) -I- A^yU^^^ -I-... 


( 6 ) 


The values of {k = 0,1,2,...) can be found by an iterative Newton-Raphson optimiza¬ 
tion of the occupation error with respect to the chemical potential using the relation 


1 dP \ 

Xk Q^ik) J 


= Py, = y5P(°)(/-P(°)), 

A=0 


(7) 


which for the approximate expanded density matrix, Eqs. ([^ and (|^, is exact in the limit 
M —)■ oo. The trace of P^, dehned here, gives the change in occupation with respect to a 
change in yU. The small deviation from the exact analytic derivative for a hnite expansion 
order M is in practice insignihcant, though for very low values of M the rate of convergence 
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will be slightly lower than quadratic in analogy to quasi Newton schemes. In combination 
with low temperatures, low values of M may also lead to loss of convergence (see Table 
0 However, in this case we could typically use regular zero temperature response theory, 
or alternatively, a modihed search routine to adjust for the correct occupation would be 
needed. 

The canonical density matrix perturbation theory based on Eqs. ([T]j^ above, which is 
our Erst key result, is summarized by Algorithm [T] for up to third order response. Each 
inner loop requires the solution of a system of linear equations, which can be achieved with 
the conjugate gradient method using as initial guesses. The linear conjugate gradient 
method |19] is ideal for this purpose, since it efficiently can take advantage of matrix sparsity 
to reduce the scaling of the computational cost na- Generalizations and modihcations to 
higher order response, grand canonical schemes (with a hxed value of /i), or spin-polarized 
(unrestricted) systems are straightforward. It is interesting to note that the system matrices 
on the left hand-side of the inner loop of Algorithmare all the same, i. e. The same 

inverse of would therefore give the response Xn^'^ for all orders k. The conditioning of 
the response algorithm should therefore be the same as for the original Oth-order expansion. 
The system matrix is very well conditioned with a spectral condition number smaller 
than or equal to 2 |18] at any point of the algorithm. In the limit of low temperature and 
high n, —)■ / and in the opposite limit of high temperatures does the condition number 

go to 1 as 1/2. The well behaved conditioning is independent of the condition 

number of the Hamiltonian used in the initialization. 

III. FREE ENERGY RESPONSE THEORY 

To study the quantum response valid at hnite electronic temperatures, the electronic 
entropy contribution to the free energy has to be considered. We will look at two different 
situations: a) non self-consistent band energy response as in regular tight-binding theory 
using an orthogonal matrix representation and b) self-consistent free energy response as in 
density functional or Hartree-Fock theory using a non-orthogonal formulation. To clearly 
separate the two cases we will use two different notations. For the orthogonal tight-binding 
like formulation we keep using H and P, which is consistent with the previous sections, and 
for the self-consistent free energy response we use F and D for the non-orthogonal matrix 
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Algorithm 1 Canonical density matrix response theory 
M ^ Number of recursion steps 

^(0) ^ Initial guess 

ph) ^ 0 Initial guess {i = 1, 2, 3} 

(3 = l/(/cBTe) ^ Choose temperature 

while Occupation error > Tolerance do 

= ^I- 2-(2+^)/3(i/(o) - ^(0)/) 

{i = 2,3} 

for n = 1 : M do 


solve for 

-^n ? 


= 0, 

1,2,3} 

7^(0) v(0) 

_ /^(O) 

— '-'n-l 





7^(0) v(l) 

_ ^(1) 

+ 

^(1) 

^n-1 

y(0) 


7^(0) v(2) 

II 

+ 

^(2) 

^n-1 

y(0) 

-^n 


7^(0) v(3) 

II 

+ 

^(3) 

^n-1 

y(0) 


end for 






II 

{i = 0 

,1: 

,2,3} 




^(0) = //(o) + (iVe - Tr[p(°)])/Tr[P^] 

_ rr[P«]/Tr[P^], {i = 1, 2, 3} 

Occupation error = |Tr[p(°)] — Ne\ + |T'r[ph)]| 

end while 
using: 

= /3pO)(/-p(o)) 

= 2x^^\xi^'> -1) +1 

cP {i.i>0, m = 0.1,2,3} 

B!r>=2(xP-CP), {m = 0,1,2,3} _ 

representations and F-^ and D-^ for the orthogonalized representations, as is explained in 
the sections below. 
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A. Non self-consistent tight-binding-like free energy response 


In a simple tight-binding like formulation, the expansion terms for the canonical free 
energy, 


ff(A) = Tr[PiX)HiX)] - T,5[P(A)] = 

= ffW + + A2fi(2) + . . _ ^ 

generated by a perturbation in H{X), Eq. Q, with the electronic entropy 


( 8 ) 


are given by 


S[P] = -kBTr[Pln{P) + (/ - P) ln(J - P)], 


Q{m) ^ 


(9) 


( 10 ) 


k=l 


This expression, with P^^'^ calculated from our canonical density matrix perturbation scheme 
in Algorithm is a straightforward generalization of the conventional Tg = 0 limit of 
the “n + 1” rule [19] and follows directly from the fact that the hrst order response term 
j^^[^(o)p{i)] jg cancelled by the response in the entropy |16]. Higher-order derivatives of 
order n + 1 therefore contain at most a derivative of order n in the density matrix. This 
generalization is possible only by including the entropy term in Eq. ([^, which is required to 
provide a variationally correct description of the energetics. We have not been able to hnd 
any explicit density matrix expressions for Wigner’s 2n +1 rule [THl UHl EU-ES] that are valid 


also at hnite temperatures. A more detailed derivation of Eq (10) is given in the appendix. 


B. Self-consistent free energy response 

In self-consistent hrst principles approaches such as Hartree-Fock theory [SSI (density 
functional and self-consistent tight-binding theory, although different, follow equivalently) 
the free energy in the restricted case (without spin polarization) is given by a constrained 
minimization of the functional 

Hscf[P] = 2Tr[hD] + Tr[DG{D)] - 2T,S[D^], (11) 

under the condition that 2Tr[PS'] = Ne, where Ne is the number of electrons (two in 
each occupied state). Here is the orthogonalized representation of the Hartree-Fock 
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density matrix D such that D = ZD^Z^, and the orthogonalized effective single-particle 
Hamiltonian is given by F-^ = Z'^FZ, where the Fockian F = h + G{D) and Z is the inverse 
factor of the basis set overlap matrix S such that Z'^SZ = L The density matrix, D, is 
thus given by 


D 




-1 


( 12 ) 


which can be calculated through the recursive Fermi operator expansion in Eqs. ([^ and 
(|^. Here h is the usual one-electron term and G{D) is the conventional two-electron part 
including the Coulomb J and exchange term K, respectively [SB]. In density functional 
theory, the Fockian F is replaced by the corresponding Kohn-Sham Hamiltonian, where the 
exchange term K is substituted with the exchange-correlation potential term. Notice that to 
make a clear distinction to the non-self-consistent response we use the notation D and F for 
the self-consistent Hartree-Fock density matrix and Fockian, i.e. the effective single-particle 
Hamiltonian. 

With a basis-set independent first order perturbation in the one-electron term. 


h(\) = 


( 13 ) 


for example due to an external electric field, the self-consistent response in the density matrix 
is given by the solution of the coupled perturbed self-consistent field (SCF) equations as in 
density functional perturbation theory: 

F(A) = + AF)(i) + ...), 


F^(A) = Z^F{X)Z, 


(14) 


D{\) 


Z 


g/3(F-L(A)-M/) 



where D and F are expanded in terms of A, i.e. 


D{X) = + XD^^'> + A^D^^) + ... , 

F(A) = + AF(i) + A^F^^) + ... . 


(15) 


The coupled response equations above are solved in each iteration using the canonical density 
matrix perturbation theory as implemented in Algorithm with H and P replaced by F-^ 






and . At self consistency, the free energy expansion terms, 


hlsCF(A) = hlscF[-D^°^] + AhlgcF + A^hlgCF + • • • ) 


( 16 ) 


are given by 


VL 


(m) 

SCF 


= —m > 0. 


m 


(17) 


This simple and convenient expression for the basis-set independent free energy response. 


which follows (see Appendix) from Eq. (10), is another key resnlt of this paper. The free 
energy response theory presented here provides a general techniqne to perform rednced 
complexity calcnlations of, for example, temperatnre-dependent static polarizabilities and 
hyperpolarizabilities mim]. 


IV. FINITE DIFFERENCE APPROXIMATIONS 


An alternative to the canonical density matrix pertnrbation theory is to perform calcn¬ 
lations with hnite pertnrbations and nse hnite difference approximations of the free energy 
derivatives. However, this can be far from trivial becanse the nnmerical errors are sometimes 
difhcnlt to estimate and control, in particnlar for high temperatnre hyperpolarizabilities. 
Nevertheless, by nsing hnite steps 5\ of the pertnrbations in h, combined with mnlti-point 
high-order hnite diherence schemes, it is sometimes possible to reach good accnracy. This 
is illnstrated in Fig. [T| which shows the hnite diherence error in the approximation of the 
second order free energy response, flgcp, with respect to an external electric held for a self- 
consistent tight-binding model [FTHUO] as implemented in the electronic strnctnre program 
package LATTE [2HIE9]. Finite diherence calcnlations of higher-order hyperpolarizabilities 
show similar behavior. 

In a hnite diherence approximation it is difhcnlt to know a priori what step size AA to 


nse for the pertnrbations in Eq. (13). Errors may be large unless careful numerical 


testing is performed. This can be expensive and even when an optimal step size has been 
found, the computational cost is still higher than the analytical approach. For example, 
to calculate the second order response using the hve point hnite diherence scheme has a 
computational cost of about 5 times a ground state calculation, whereas the cost for the 
density matrix perturbation theory is only about 3 times larger. This cost estimate does 
not include the additional entropy calculations. The calculation of the entropy is difficult 
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Log^CA^i) 


( 2 ) 

FIG. 1: The relative error compared to the “exact” derivative in Eq. (17) for using 5 

and 9 point central difference schemes for the calculation of the second order response in the free 
energy with respect to an external electric field, i.e. the polarizability. Either the exact entropy 
expression was used, Eq. (§, or the highest order (m = 4) approximation in Eq. ( [T^ . The 
electronic temperature Tg is about 37,000 K. 


(or impossible) to perform accurately within linear scaling complexity. Computationally 
favorable formulations that are based on approximate expansions of S[P] in Eq. ([^ are 
typically poor. For example, when any of the approximate entropy expressions. 




( 18 ) 


2=1 


with the coefficients Cj(m) in Tab. [^are used, the relative error of the polarizability in Fig. 

is increased by over 6 orders of magnitude for the most accurate 9 point finite difference 
approximation. The accuracy is at best only about 0.5 percent with any of the entropy 
approximations in Eq. (18) and Tab. Only by avoiding explicit entropy calculations it 
is possible to reach a meaningful accuracy. This is possible in a finite difference approx¬ 
imation by using the finite differences of the dipole moments instead of the free energies. 
Such calculations (not shown) avoid calculating the explicit entropy term and the numerical 
accuracy is similar to the finite difference approximations using the free energies with the 
exact entropy expression as illustrated in Fig. 
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TABLE I: Coefficients for the approximate entropy expression in Eq. (18). The coefficients are 


determined from the ansatz in Eq. (18) with the requirement that the function value and a few of 


its derivatives are correct at the midpoint 0.5 of the interval [0,1] in which P has its eigenvalues. 


Ci(m) 

m = 1 

m = 2 

m = 4 

ci(m) 

41n(2) 

81n(2) - 2 

161n(2) -34/5 

C2(m) 


161n(2) - 8 

961n(2) -844/15 

C3(m) 



2561n(2) - 2336/15 

C4(m) 



2561n(2) - 2368/15 


V. FIRST PRINCIPLES RESULTS 

A. Polarizabilities and hyperpolarizabilities 

Figurej^shows the calculated temperature-dependent response for a single water molecule 
with respect to static electric helds. The calculations were performed with Hartree-Fock 
theory using the ERGO program package [361 l37] . At lower temperatures the response 
values correspond to the isotropic polarizability and hyperpolarizabilities if the values are 
multiplied by ml, i. e. the factorial of the response order. At higher temperatures this 
interpretation is less accurate because of the limited basis set description of the thermally 
excited states. For relevant temperatures below 10,000 K our calculations show a very 
small temperature dependence, which is consistent with a fairly large HOMO-LUMO gap. 
For higher temperatures the errors may be signihcant, since the Gaussian basis set used 
here (cc-pVDZ) was not designed for high-temperature expansions. The calculations were 
performed for a single molecule in the gas phase. For periodic boundary conditions the 
position and the dipole moment operator are not well dehned. In this case the techniques 
developed within the modern theory of polarizability can be applied [HTHM] . 

The response properties converges quickly as a function of the number of recursion steps 
(M) in the canonical density matrix response expansion in Alg. which is illustrated in Tab. 
[ITj At higher temperatures we see a slightly slower convergence, and at low temperatures 
and with a small number of recursion steps there can be problems with convergence of the 
occupation, since the chemical potential derivative estimate is less 
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(‘ 2 \ 1 

FIG. 2: (color online) The temperature dependent isotropic second-order response fIgQp[ 3 (xx -|- 
yy + zz)] = 3 (IIg()p[xx] -|- ng( 3 p[yy] -|- flg()p[z 2 :]), and the third order and fourth order response 
in the x direction. The canonical density matrix response Algorithm for restricted Hartree- 
Fock theory (RHF) with a Gaussian basis set (cc-pVDZ) was used. The xyz-coordinates of the 
molecule: {O (0.0,0.0, 0.0); H (-1.809,0.0,0.0); H (0.453549,1.751221,0.0)} in atomic units. 
As a comparison and validation we show five-point finite difference calculations of the free energy 
derivatives. At low electronic temperatures the second-order response corresponds to 1/2 times 
the isotropic polarizability (see Tab. 0. 

accurate. In this case we may prefer to use a regular zero-temperature response calculation. 

B. Linear scaling complexity 

It is easy to understand the potential for a linear scaling implementation of canonical 
density matrix perturbation theory. Owing to nearsightedness p dg sal m], both the 
Hamiltonian and its perturbations, as well as the density matrix and its response, have 
sparse matrix representations for non-metallic materials when local basis set representations 
are used. The number of signihcant matrix elements above some small numerical threshold 
(or machine precision) then scales only linearly with the number of atoms for sufficiently 
large systems. In this case, since all operations in the canonical density matrix perturbation 
scheme in Algorithm [T] are based on matrix-matrix operations, the computational cost scales 
only linearly with system size if sparse matrix algebra is used in the calculations. This is not 
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1 f2) 

TABLE II: Convergence of the isotropic polarizability Oiso = (^SCf[^^] + ^scpby] + 

^SCpC^y]) three different electronic temperatures Tg (1000 K, 30,000 K and 100,000 K) 
as a function of the number of recursion steps (M) in the canonical density matrix re¬ 
sponse expansion in Algorithm for a water molecule calculated from restricted Hartree-Fock 
theory (RHF) with a Gaussian basis set (cc-pVDZ). The xyz-coordinates of the molecule: 
{O (0.0, 0.0, 0.0); H (—1.809, 0.0,0.0); H (0.453549,1.751221, 0.0)} in atomic units. As a compar¬ 
ison and validation we show the Te = 0 K result of the isotropic polarizability, which was calculated 
by solving the linear response time-dependent Hartree-Fock (or RPA) equations [35] as implemeted 
in the ERGO program package [36l EZ] applied for the zero-frequency case. 


Te (K) 

M 

Q^iso 

Te (K) 

M 

*^iso (S'*!!.) 

Te{K) 

M 

Q^iso 

0 (ERGO) 

n/a 

-5.0112528623 







1000 

6 

no convergence 

40,000 

6 

-6.8540449154 

100,000 

6 

-7.5204026148 

1000 

8 

-5.0112527697 

40,000 

8 

-6.8538983381 

100,000 

8 

-7.5198385798 

1000 

10 

-5.0112527697 

40,000 

10 

-6.8538891617 

100,000 

10 

-7.5198033131 

1000 

12 

-5.0112527697 

40,000 

12 

-6.8538885881 

100,000 

12 

-7.5198011089 

1000 

14 

-5.0112527697 

40,000 

14 

-6.8538885522 

100,000 

14 

-7.5198009711 

1000 

16 

-5.0112527697 

40,000 

16 

-6.8538885500 

100,000 

16 

-7.5198009625 


possible in regular Rayleigh-Schrodinger perturbation theory, which requires the calculation 
of individual eigenvalues and eigenfunctions. Figurej^shows the number of non-zero elements 
above threshold as a function of system size for the density matrix and its first and second 
order response with respect to an electric dipole perturbation. The test systems are simple 
one-dimensional hydrocarbon chains of various lengths and the calculations where performed 
based on Hartree-Fock theory using a small Gaussian (STO-3G) basis. Gaussian basis 
sets were not designed for the high-temperature expansions demonstrated here and we can 
expect that the accuracy is limited. The simulations therefore only serve as a schematic 
demonstration of the expected behavior. For example, at higher temperatures the locality, 
i.e. the matrix sparsity, is increased similar to what is found for larger HOMO-LUMO 
gaps and for higher order response the locality decreases, as has been seen in previous 
studies of the zero temperature case miiiB]. Using a larger Gaussian basis set should not 
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FIG. 3: (color online) The sparsity scaling as a function of system size of the density matrix 
and its first and second order response with respect to an electric (static) dipole perturbation for 
two different electronic temperatures. The graphs show the number of non-zero elements of the 
orthogonal density matrix after a numerical threshold of 10“^. 

change this general behavior of the locality and the resnlts wonld still be nncertain. Matrix 
sparsity may also snffer and nnmerical problems may arise dne to ill-conditioning from linear 
dependencies between many Ganssians. However, this will not affect the conditioning of the 
canonical density matrix response algorithm, Alg. [T| and the low spectral condition nnmber 
of which is always < 2, bnt it wonld affect the congrnence transformation from the 

non-orthogonal atomic orbital representation of F to F-^. The inpnt data of the response 
algorithm wonld thns be less accnrate. Localized nnmerical atomic orbital basis sets that 
have been tailored specihcally for high-temperatnre expansions (and with low condition 
nnmbers of the overlap matrix) wonld then be a more appropriate choice. 


VI. SUMMARY 

In snmmary, we have presented a canonical single-particle density matrix pertnrbation 
scheme that enables the calcnlation of temperatnre dependent qnantnm response properties. 
Since onr approach avoids the calcnlation of individnal eigenvalnes and eigenfnnctions as 
well as the entropy, the theory is well adapted for rednced complexity calcnlations with a 
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computational effort that scales only linearly with the system size. However, we may expect 
very fast parallel performance also for smaller systems in the limit of dense matrix algebra, 
since the computational kernel is centered around matrix-matrix multiplications that often 
can reach close to peak performance on modern hardware. The perturbation scheme should 
be applicable to a number of existing program packages for linear scaling electronic structure 
calculations. 
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VIII. APPENDIX 


A. Recursive Fermi operator expansion 

There are several techniques to calculate matrix exponentials. For example, if we start 
with 


= 

a hrst order Taylor expansion gives 

= lim 


ga:/(2n) 

g-x/(2n) 


2n + X 


(19) 


( 20 ) 


n-5>oo \2n — X ^ 

Using this expansion we can approximate the Fermi-Dirac distribution function, <h((r), with 

(2n — xY 


<h(a;) = -1- 1] ^ = lim 


n^co (2n -|- xY + (2n — xY ’ 


such that 


<h(2n — Anx) 


X 


x^ + (1 -xY' 


( 21 ) 


( 22 ) 
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which is accurate for large values of n. The Fade polynomial function 


fn{x) ~ n i \n 

+ [1 — xp 

can be expanded recursively, since 

fmxn{x) /m(/n(^))- 


( 23 ) 


(24) 


This particular property enables a rapid high-order expansions in only a few iterations in 
the recursive Fermi-operator expansion, 


$ [I3{ei - /i)] = <F(2n - AnXi) 


(25) 


~ ■ ■ f2{Xi) . . .)), 

where 

( 26 ) 

with the recursion repeated m times, i. e. for n = 2™'. In 30 steps (m = 30) this gives an 
expansion order of the Fade polynomial of over 1 billion, but often less than 10 steps are 
needed. 

The density matrix at finite electronic temperatures, 

P ^ + xj ^ (27) 


can now be calculated with the recursive grand canonical Fermi operator expansion, 

P = /2 (. .. /2 (/2 (^1/ - 2-<=+“p(/f - f,/) j (28) 

which forms the starting point in Eq. ([^, with J^niX) = / 2 (X). The recursive grand 
canonical Fermi operator expansion, derivations, convergence analysis, and tests with various 
basis sets have been published previously in Refs. 


B. Perturbation response for the non-self-consistent single particle free energy 


To derive Eq. (10) we start by noting that from the definition of the density matrix 
response and the perturbations in the Hamiltonian, Eqs. (|^ and (|^, we have 

Qk 




P(A) 


_ p[k] _ }^\p{k) 


(29) 


A=0 


16 






and 


Qk 

dX^ 


H{X) 


= 


(30) 


A=0 


where we use square brackets for the regular Taylor expansion terms, and and 
round brackets, and for the perturbation expansions as in Eqs. (|^ and (|^. 

Thereafter we can calculate the response terms from the derivatives of the free energy 
expression in Eq. ([^, i. e. 

^ fi'^ 

(31) 


1 


ml dX"^ 


A=0 


It is easy to see that the first derivative of the entropy term iS[P(A)] in Eq. ([^ is given by 
lS[P{X)]\^^^ = -ksTr [(ln(P) - ln(/ - P))PW] 


= -keTr [in (P(J - P)-i)pW] = -ksTr [ln(e-^('f^-^'f))pW] 
= fcB/3Tr[PPb)], 


(32) 


since we have a canonical perturbation Tr[p[^]] = 0 and P = + /] ^. This means 

that the hrst order response in the free energy hi (A) is given by 


fi(i) = |-(](A)|^^^ = Tr[p(i)pW] + Tr[p(°)p(i)] - Tefcs/3Tr[PP(i)] 
= Tr[pWp[°l] = Tr[P(^)p(°)]. 

For the second order expansion we find that 

= lTr[pWp[il + p-[2]p[o]] 

= I (Tr[P(i)p(i)] + 2Tr[P(2)p(o)]) _ 

For the third order expansion we find that 

fi(3) = + ffPIplil + J/Plplol + 

= |rr[2J?<‘>FPI + + 62?P>P<“> + 22?<2lFl'l] 

= i (rr|ffl‘)pP>] + 2rr[J/(^lF(‘)] + 3rr|/fl=»Pl“>]) . 


(33) 


(34) 


(36) 


The straightforward mth-order generalization from consecutive derivatives gives Eq. (10). 


C. Basis-set independent self-consistent free energy response 


To derive the basis-set independent response of the free energy in Eq. (17) we first cal¬ 
culate the hrst order derivative of 


nscF[D] = 2Tr[hD] + Tr[PG(P)] - 2T,5[P^], 


(36) 
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with respect to A in Eq. (13), i. e. 


= 2Tr[h(i)E) + + 2Tr[D^^^G{D)] - 2TefcB/STr[F^E)^'^'] 

= 2Tr[h^^^D] + 2Tr[{h + G(E)))E)W] - 2Tr[F^D^^^\ (37) 

= 2Tr[h(i)E)[o]] + 2Tr[FE)W] - 2Tr[Z^ 

= 2Tr[h(i)E)[°]] + 2Tr[FE)W] - 2Tr[FE)W] = 2Tr[h(^)E)M] 


where we have derived the entropy derivative as in Eq. (32) above, used the definition of the 
Fockian, F = h + G{D), applied the congruence transformation between the orthogonal and 
non-orthogonal representations, e. g. F^ = Z'^FZ, and the cyclic permutation under the 
trace. Subsequent derivatives, analogous to the previous Appendix subsection above, gives 
Eq. (|g. 
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